home *** CD-ROM | disk | FTP | other *** search
/ CD Concept 6 / CD Concept 06.iso / mac / UTILITAIRE / RLaB / testmatrix / augment.r next >
Text File  |  1994-12-20  |  2KB  |  63 lines
  1. //-------------------------------------------------------------------//
  2.  
  3. // Synopsis:    Augmented system matrix.
  4.  
  5. // Syntax:    augment ( A )
  6. //        augment ( A , ALPHA )
  7.  
  8. // Description:
  9.  
  10. //    augment(A) is the square matrix [EYE(m) A; A' ZEROS(n)] of
  11. //    dimension m+n, where A is m-by-n.  It is the symmetric and
  12. //    indefinite coefficient matrix of the augmented system
  13. //    associated with a least squares problem minimize NORM(A*x-b).
  14.  
  15. //      Special case: if A is a scalar, n say, then AUGMENT(A) is the
  16. //    same as AUGMENT(RANDN(p,q)) where n = p+q and p = ROUND(n/2),
  17. //    that is, a random augmented matrix of dimension n is produced. 
  18.  
  19. //    If a second, scalar argument ALPHA is supplied, then the (1,1)
  20. //    block of A is ALPHA*EYE(m). The eigenvalues of AUGMENT(A) are
  21. //    given in terms of the singular values s(i) of A (where m>n) by
  22. //    1/2 +/- SQRT( s(i)^2 + 1/4 ),  i=1:n  (2n eigenvalues), 1,
  23. //    (m-n eigenvalues). 
  24.  
  25. //      If m < n then the first expression provides 2m eigenvalues and
  26. //    the remaining n-m eigenvalues are zero.
  27.  
  28. //      See also SPAUGMENT.
  29.  
  30. //      Reference:
  31. //      G.H. Golub and C.F. Van Loan, Matrix Computations, Second
  32. //      Edition, Johns Hopkins University Press, Baltimore, Maryland,
  33. //      1989, sec. 5.6.4.
  34.  
  35. //    This file is a translation of augment.m from version 2.0 of
  36. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  37. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  38.  
  39. //-------------------------------------------------------------------//
  40.  
  41. augment = function ( A, alpha )
  42. {
  43.   local (A, alpha)
  44.  
  45.   m = A.nr; n = A.nc;
  46.  
  47.   if (!exist (alpha)) { alpha = 1; }
  48.  
  49.   if (max(m,n) == 1)
  50.   {
  51.     n = A;
  52.     p = round(n/2);
  53.     q = n - p;
  54.     rand("normal", 0, 1);
  55.     A = rand(p,q);
  56.     m = p; 
  57.     n = q;
  58.   }
  59.  
  60.   C = [alpha*eye(m,m), A; A', zeros(n,n)];
  61.   return C;
  62. };
  63.